
% Driver script for the function computererrbounds_network()
% 
% P. Ajith, 8 Aug 2009
% 
% $Id:$

clear 

tic 

Params.startPhase = 0;
Params.startTime = 0;
Params.fLower = 10;
Params.eta = 0.25;
Params.totalMass = 20;

detNameVec = {'LHO','LLO','Virgo'};
noisePSDVec = {'AdvLIGO','AdvLIGO','AdvVirgo'};
%lambdaVec = {'dL','M','eta','theta','phi','t0','phi0','psi'};
lambdaVec = {'M','eta','theta','phi','dL','t0','phi0','psi'};
projectOutDims = 1;

Params.theta = 0.810795263059348;
Params.phi = -2.08405676907401;
Params.psi = 0;
Params.iota = 0;
Params.dL = 1000;

% grid points of theta and phi
nGrid = 30;
thetaVec = linspace(0,pi,nGrid);
phiVec = linspace(-pi, pi, nGrid);

% index of the sky-location parameters (theta, phi) in the lambdaVec
thetaIdx = find(strcmp(lambdaVec,'theta'));
phiIdx = find(strcmp(lambdaVec,'phi'));

% compute the error bounds for each pair of masses
for iTheta=1:length(thetaVec)
        
    Params.theta = thetaVec(iTheta);
    for iPhi =1:length(phiVec)

        fprintf('iPhi = %d iTheta = %d\n', iPhi, iTheta);

        Params.phi = phiVec(iPhi);

        % compute the error bounds at this sky location using a network 
        % of detectors
        [sigmaVec, covMat, snr(iTheta, iPhi), snrVec, condNum(iTheta, iPhi), condNumVec...
            ] = computeerrbounds_network(Params, detNameVec, noisePSDVec, lambdaVec, ...
            projectOutDims);

        % save the RMS error in each parameter in a matrix in a structure
        % Param(iParam).sigmaMat
        for iParam = 1:length(sigmaVec)
            Param(iParam).sigmaMat(iTheta, iPhi) = sigmaVec(iParam);
        end

        % extract the covariance-matrix element describing the
        % correlation btween theta and phi
        corrPhiTheta(iTheta, iPhi) = covMat(phiIdx,thetaIdx);

        % also keep a matrix version of the sky-position angles
        phiMat(iTheta, iPhi) = phiVec(iPhi);
        thetaMat(iTheta, iPhi) = thetaVec(iTheta);
    
    end
end

networkName = strcat(detNameVec{1:end});

% plot the network SNR
fH = figure;
imagesc(phiVec*180/pi,thetaVec*180/pi, snr)
colorbar
colormap('hot')
grid on
xlabel('\Phi (deg)')
ylabel('\Theta (deg)')
title(sprintf('Network SNR (%s)', networkName))

% save figure
figName = sprintf('SNRSkyMap_%s_dL%dMpc_nGrid_%d', networkName, Params.dL, nGrid);
saveas(fH, figName, 'fig')

% compute the RMS error in dOmega
sigmaTheta = real(Param(thetaIdx).sigmaMat);
sigmaPhi = real(Param(phiIdx).sigmaMat);
dOmega = 2*pi*sin(thetaMat).*sqrt(sigmaTheta.^2.*sigmaPhi.^2 - corrPhiTheta.^2);
dOmega = dOmega*(180/pi)^2;     % in square degrees

dataParamVec = lambdaVec;
dataParamVec{end} = 'Omega';
Param(end+1).sigmaMat = dOmega;
Param(1).snr = snr;
Param(1).lambdaVec = dataParamVec;
Param(1).corrPhiTheta = corrPhiTheta;
Param(1).phiMat = phiMat;
Param(1).thetaMat = thetaMat;

fH = figure;
imagesc(phiVec*180/pi,thetaVec*180/pi, log10(real(dOmega)))
colorbar
colormap('hot')
grid on
xlabel('\Phi (deg)')
ylabel('\Theta (deg)')
title(sprintf('log \\Delta \\Omega [sq. degree] (%s)', networkName))
    
% save figures and the data as mat files
figName = sprintf('ErrSkyMap_dOmega_%s_dL%dMpc_nGrid%d', networkName, Params.dL,nGrid);
saveas(fH, figName, 'fig')
save(figName, 'Param')

% plot the RMS errors of all parameters
for iParam = 1:length(sigmaVec)

    fH = figure(iParam+2);
    imagesc(phiVec*180/pi,thetaVec*180/pi, log10(real(Param(iParam).sigmaMat)))
    colorbar
    grid on
    colormap('hot')
    xlabel('\Phi (deg)')
    ylabel('\Theta (deg)')
    title(sprintf('log \\Delta %s (%s)', lambdaVec{iParam},networkName))

    % save figures
    figName = sprintf('ErrSkyMap_%s_%s_dL%dMpc_nGrid%d', ...
        lambdaVec{iParam}, networkName, Params.dL, nGrid);
    saveas(fH, figName, 'fig')
end

toc
